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Abstract 


The model Tip4p/e for water is tested for the presence of thermodynamic and dy¬ 
namic anomalies. Molecular dynamic simulations for this model were performed and 
we show that for this system the density versus temperature at constant pressure ex¬ 
hibits a maximum. In addition we also show that the diffusion coefficient versus density 
at constant temperature has a maximum and a minimum. The anomalous behavior 
of the density and of the diffusion coefficient obey the water hierachy. The results for 
the Tip4p-£ are consistent with experiments and when compared with the Tip4p-2005 
model show similar results a variety of physical properties and better performance for 
the dielectric constant. 


Introduction 

Water is a fascinating molecule. Even though present in our everyday life, it shows a number 
of properties that are still not well described. 1 2 For example, most liquids contract upon 
cooling. This is not the case of water, a liquid where the specific volume at ambient pressure 
starts to increase when cooled below 4 C at atmospheric pressure.^ 3 In addition, in a certain 
range of pressures, water also exhibits an anomalous increase of compressibility and of the 
specific heat upon cooling. 4 6 Water also has dynamic anomalies. Experiments show that 
the diffusion constant, Z), increases on compression at low temperature, T, up to a maximum 
D ni ax(T) at p — pomax{T). The behavior of normal liquids, with D decreasing on compression, 
is restored in water only at high p , e.g. for p > pDmax ~ 11 kbar at 10°C . 7 

In addition to the measured anomalies of water, theoretical analysis predicted anoma¬ 
lies 2 8 that are located in regions of the pressure versus temperature phase diagram of difficult 
access experimentally. Consequently, simulations became an interesting tool to test these 
theories. Then the challenge faced when developing a computational strategy is to design a 
model that would be general enough to describe the different behaviors of water and simple 
enough to be computationally treatable. The later prerequisite at the moment preclude the 
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consideration of quantum effects and polarization. Both polarization and quantum effects, 
however, seems to play a relevant role in the anomalous properties of water particularly 
when charges and interfaces are present. In order to circumvent this difficulty without loos¬ 
ing the simplicity required for simulation purposes, a number of atomistic models has been 
developed with the assumption that polarization and quantum effects were included in an 
averaged way. 

These atomistic models are characterized by representing the charges in water by two, 
three, four or even five points. Then, the interactions were modelled by a classical Lennard- 
Jones for the hardcore interactions and the electrostatic interactions for the charges. This 
leads to the following parameters that need to be specified: the values and positions of the 
charges, the positions and masses of the atoms and the energy and size for the LJ interaction. 
Then, the crucial step in the modelling process is the choice of the set of quantities used to 
fit these parameters. This set should be small but appropriated to guarantee that the model 
reproduces as many properties of water as possible at least in a certain range of temperatures 
and pressures. 

Within the non-polarizable models the 4-site form represented an advance. It was first 
proposed by Bernal and Fowler 9 along with a set of parameters based on calculations for 
properties of the monomer, dimer, and ice. The fours points are the position of the oxygen 
and hydrogens and the location , M , of the negative charge. Within the Tip4p each hydrogen 
carries a positive charge, qn, while the negative charge, qM is located at a position vom from 
the oxygen between the two hydrogens. The angle between the oxygens and the hydrogens, 
104.52°, and the distance between the oxygen and the hydrogen, toh = 0.9572, where fixed 
to reproduce the ice structure. The Bernal and Fowler model, however, gives a poor results 
for the liquid properties at room 25°C and atmospheric pressure. A reparametrization of 
this model gave rise to the TIP4P 10 model that shows good agreement with the density 
at 25°C and 1 atm and an excellent value for the vaporization enthalpy. In addition this 
model provides a reasonable description of some solid phases and reproduces qualitatively 
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the phase diagram 11 14 while the results for the SPC/E and TIP5P models are quite poor. 
Unfortunately it gives a value that is too low for the temperature of maximum density and 
of melting. Then, it became clear that good model of water should provide the behavior 
of the liquid, particularly the value of the density anomaly at atmospheric pressure, and a 
reasonable description of the solid phases. For that purpose the TIP4P/2005 12 was created. 

It was designed to match the density at the temperature of maximum density but yields a 
slightly low melting temperature and a somewhat large vaporization enthalpy. 

In order to test the TIP4P/2005 model against other options, Vega et al. 14 have compared 
a number of the non-polarizable models. The strategy was to select a set of water properties 
and compared the results obtained by different models with the experiments. They estab¬ 
lished that the best model to reproduce the properties they have selected is the four sites 
TIP4P/2005 12 13 followed by the three sites SPC/E model. 15 The only drawback of these 
models is that they do not give a good description for the dielectric constant of water. Since 
the dielectric constant is fundamental for understanding the behavior of mixtures of water 
and other substances, particularly polar molecules, these water models are not appropriated 
to analyse these mixtures. 

In order to circumvent this difficulty without loosing the advantages of the TIP4P/2005,^D31 
Fuentes et al. 17 * developed the non-polarizable TIP4P/e rigid model. This potential is 
parametrized to give the experimental value of the density and of the dielectric constant 
at 4 °C and atmospheric pressure. This new model showed that it is capable of reproducing 
some thermodynamic quantities* 77 * obtained by the TIP4P/2005. 12 13 In addition it gives a 
good agreement with the experiments for the isothermal compressibility and dielectric con¬ 
stant at different pressures and temperatures what is not observed in the non-polarizable 
models. 

In addition to the thermodynamic anomalies, water also show a singular mobility. While 
experimental results show that the diffusion coefficient of water decreases with decreasing 
pressures up to crystallization, simulations with SPC/E water show that this system if kept 
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liquid reaches a minimum 18 "2S1 a t negative pressures. Then the pressure and the temperature 
of the maximum and the minimum of the mobility define a region of diffusion anomaly. 
This region englobes the density anomaly defining the hierarchy of the anomalies.® 20 This 
hierarchy has been employed to conceptualize the mechanism behind the thermodynamic 
and dynamic unusual behavior of water. This result suggests that the thermodynamic and 
the dynamic anomalies are not independent, but are related by the competition of two length 
scales: bonding and non bonding.^ 

Therefore it would be desirable that a model for water would be capable to capture not 
only the thermodynamic anomalies but also the dynamical anomalous behavior of water. In 
this paper we test if the TIP4P/C model also shows the dynamic anomalous region in the 
pressure versus temperature phase diagram observed in water. We compute the diffusion 
coefficient, D , versus temperature for various densities and temperatures. Then the location 
in the pressure versus temperature of the maximum and minimum of the diffusion coefficient 
D are compared with the density extrema and checked if the TIP4P/e has the hierarchy 
observed in the experiments. Our results are also compared with the TIP4P/2005 12 14 
model. Finally a summary of the thermodynamic, dynamic and structural properties of this 
model is compared with experiments and with the TIP4P/2005 mode in the spirit of the 
grading proposed by Vega et ah 14 

The remaining of this paper goes as follows. In the section 2 the force fields for the 
TIP4P/2005 and TIP4P/e are presented. In the section 3 the simulations are explained and 
in section 4 results are analysed. Conclusions are shown in section 5. 

The Models 

The different propositions for the four site models have in common the format illustrated in 
the Figure [lj The system is represented by the oxygen and the two hydrogens. The distance 
between the oxygen and the hydrogen is given by ron = 0.9572 while the angle between the 
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Figure 1: Model of four points for water: one oxygen, two hydrogens and a point M where 
the negative charge is located. 

two hydrogens and the oxygen is Qhoh — 104.52°. These quantities were set in order to 
give the appropriated ice form. The oxygen has mass mo = 15.999 g/mol and the hydrogens 
have mass mo = 1.008 g/mol. The shared electrons between the hydrogens and the oxygen 
are closer to the oxygen. This is represented by the two positive charges, qn, one at each 
hydrogen and a negative charge qM ~ 2 qn at a point M distant roM from the oxygen and 
located between the two hydrogens. Each water molecule i has a kinetic energy, Kj given by 

Ki = ^m 0 vf (1) 


where v ; - is the velocity of the molecule. Two water molecule interact through a potential with 
two contributions, a Lennard-Jones (LJ) between the oxygens and electrostatic interactions 
between the hydrogens and the negative charge at the point M, namely 


u(rij) = 4e 0 o 


f f Goo\ 1 y-i y-i qtqj 

\ r ij ) V r ij ) \ 4n£ o h\ h\ r d 


( 2 ) 


where is the distance between atom i and j, qt is the electric charge of atom i, £o is the 
permittivity of vacuum, £qo is the LJ energy scale of the oxygen interactions and Oqo th e 
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diameter for an 00 pair. The model has one LJ site and charge on the oxygen atom and 
additionally a charge on every hydrogen atom. 

In this paper some thermodynamic, dynamic and structural properties of two force fields 
are compared: TIP4P/2005 and the TIP4P/e. For the first model the parameter, roMi &00- 
qn, (<?M = 2 qn) and Ooo are selected by imposing that the model reproduces the maximum 
density at T — 4°C and atmospheric pressure. For the TIP4P/e the parameters are chosen 
so the model not only reproduces the density but also the dielectric contant at T — 4 °C at 
atmospheric pressure. The parameters of these two force fields are given in the Table [lj 

Table 1: Force field parameters for water models. The charge in site M is qM — 
-(2 <?//)• 


Model 

qii/e 

Qm/^ — 'ZqH/e 

r OM / 

Coo/A 

(fiooAtf)/ K 

TIP4P/C 

0.9572 

1.054 

0.105 

2.4345 

93 

TIP4P/2005 

0.9572 

1.1128 

0.1546 

2.305 

93.2 


Simulation details 

All the simulations in this work have been done for a system of 500 molecules and employing 
molecular dynamic simultions in the NVT ensemble with the package GROMACS 4.5.® The 
equations of motion are solved using the leap-frog algorithm® 24 and the time step used was 
2 fs. The Lennard-Jones potential has been switched from 10 up to a cut-off distance of 10 
. Long range corrections were applied to the Lennard-Jones part of the potential (for both 
the energy and pressure). 

Ewald summations were used to deal with electrostatic contributions. The real part of 
the Coulombic potential is truncated at 10 . The Fourier component of the Ewald sums was 
evaluated by using the particle mesh Ewald (PME) method 2 ® using a grid spacing of 1.2 and 
a fourth degree polynomial for the interpolation. The simulation box is cubic throughout 
the whole simulation and the geometry of the water molecules kept constant using the shake 
procedure. 2 ' Temperature has been set to the desired value with a Nose Hoover thermostat.®^ 
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The diffusion coefficient is calculated using the mean-square displacement averaged over 


different initial times,namely 


(Ar(,) 2 } = (K'0 + <)-rM 2 }. 


( 3 ) 


From Eq. (??), the diffusion coefficient may be obtained as follows: 


D = lim (A r{t) 2 ) / 6 1 . 


( 4 ) 


The static dielectric constant is computed from the fluctuations^® of the total dipole 
moment M, 



where kg is the Boltzmann constant and T the absolute temperature. The dielectric con¬ 
stant is obtained for long simulations at constant density and temperature or at constant 
temperature and pressure. 

Results 

First, the temperatures of maximum density for the different pressures were computed for 
both TIP4P/2005 and TIP4P/e models as follows. In the NVT ensemble this is done by 
relating the minimum of the isochores at the pressure versus temperature phase diagram. 
Using the Maxwell relation, 



{ df )v( 7p )T 



( 6 ) 


the maximum of p(T) versus temperature at constant pressure given by ( dp/dT)p — 0 is 
equivalent to the minimum of ( dP/dT) p = 0. While the former is suitable for NPT-constant 
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Figure 2: Pressure versus temperature phase diagram. The red circles are the isochores for 
the TIP4P/2005 model while the blue squares are the isochores for the TIP4P/e model. The 
minimum of the isochores show the TMD location for each model. The black line indicates 
the experimental results for the TMD. indicate isochores.^ 


experiments/simulations the latter is more convenient for onr NVT-ensemble study, thus 
adopted in this work. 

Figure [2] illustrates the pressure versus temperature phase digram where the isochores 
for the TIP4P/2005 and for te TIP4P/e models are shown as red circles and blue squares 
respectively. The minimum of the isochores are also illustrated as red losangles for the 
TIP4P/2005 model and blue triangles for the TIP4P/e model. These lines locate the Tem¬ 
perature of Maximum Density (TMD). The simulations give a good agreement with the 
experimental results for the TMD represented by a black solid line. The two models are 
quite equivalent for the location of the TMD what is not surprising since both are adjusted 
to give the location of the density maximum at atmospherica pressure. 

In addition to the thermodynamic anomaly the diffusion anomaly is also analyzed. Fig- 
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Figure 3: Diffusion coefficient versus density for various temperatures. The red circles illus¬ 
trated D for the TIP4P/2005 model while the blue circles show the TMD for the TIP4P/C 
model. Both maximum of D. 


ure [3] shows the diffusion coefficient versus density for both models for a range of temper¬ 
atures. For both models de diffusion coefficient versus density graph has a maximum and 
minimum for various temperatures. The values of the temperature of the maximum and 
minimum diffusion coefficients are consistent with the values obtained by experiments . 1 The 
pressures for the maximum and minimum D, however, give higher values when compared 
with the experimental results^ what might be attributed to the rigidity of both models. 

The hierarchy of anomalies of the TIP4P/£ model shown in the Figure ?? confirms the 
predicted behavior that the region in the pressure versus temperature phase diagram in which 
the diffusion anomaly is present involves the region where the density anomaly appears.^ 

We also check the behavior of the dielectric constant, £, with the temperature and density. 
Figure [5] shows £ as a function of the temperature for different densities for the TIP4P/2005 
(red circles) and TIP4P/£ (blue squares) models. The TIP4P/2005 model shows much lower 
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Figure 4: Pressure versus temperature phase diagram illustrating the temperature of maxi¬ 
mum density (black filled circles) and the maximum (blue filled squares) and minimum (blue 
empty squares) of the diffusion coefficient for the TlP4P/e model. 

values of dielectric constant when compared with the TIP4P/£ model and experiments.® The 
fact that the TIP4P/C gives a good estimate for the dielectric constant at room pressure and 
temperature is not surprising since the model was fitted to give this result. It is interesting 
to observe that for this result is also preserved for other values of temperatures. 

Since the TIP4P/e does not provide a good evaluation of the pressure for the maximum of 
the diffusion anomaly it is important to verify if at high pressures the dielectric constant fails 
to agree with the experiments. Figure 7 illustrates the dielectric constant versus pressure 
for different temperatures for the TIP4P/e model and experiments for T — 300 K. The 
agreement between simulations and experiments is good for the low pressures. Experimental 
data for higher pressures still need to be further explored. 

In recent years Vega et al. 14 proposed to evaluate the performance of water models by 
a measure. The models received a grade from zero to ten by checking how a finite group of 
properties from the liquid, solid and gas phases predicted by the model agree with experi¬ 
mental results. In addition to equilibrium thermodynamic properties, the measure includes 
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Figure 5: Dielectric constant versus temperature for different densities from 0.75g/cm 3 to 
1 g/cm 3 (from bottom to top) for the TIP4P/2005 (red circles) and for the TIP4P/C (blue 
squares) models. The experimental 31 results (black solid line) are shown for 1 g/cm 3 . 


dynamic properties and phase transition predictions. In order to answer to the logical criti¬ 
cism that the TIP4P/e model performes well in computing the dielectric constant but might 
fail in other properties in which TIP4P/2005 gives good agreements with experiments,^ the 
measured proposed by Vega et al. 14 Iwas computed for the TIP4P/e model. 

Table ?? shows the performance of the TIP4P/e model compared with the performance 
of the TIP4P/2005 for the properties proposed by Vega et al. 14 Our results indicate that 
TIP4P/C gives good results not only for the dielectric constant, density and diffusion anoma¬ 
lies but also for the selected properties illustrated in the table. 
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Table 2 


Property 

Experiment 

data 

Quantity 

TIP4P/05 

Quantity 

TIP4P/C 

% Tol. 

Score 

TIP4P/05 

Score 

TIP4P/C 

Enthalpy of phase change / kcal mol 1 1 

AH melt 

1.44 

1.16 

1.24 

5 

6 

7 

AH vap 

10.52 

11.99 

11.74 

2.5 

4 

5 

Critical point properties I 

T c /K 

647.1 

640 

675.45 

2.5 

10 

8 

Pel g cm - 3 

0.322 

0.337 

0.2993 

2.5 

8 

7 

Pc/bar 

220.64 

146 

136 

5 

3 

2 

Surface tension/mN m 1 j 

Y500K 

71.73 

69.3 

69 

2.5 

9 

8 

Y450K 

42.88 

41.8 

43.8 

2.5 

9 

9 

Melting properties 1 

T m /K 

273.15 

252 

240 

2.5 

7 

5 

Pliq/g cm -3 

0.999 

0.993 

0.994 

0.5 

9 

9 

Psol/g cm-3 

0.917 

0.921 

0.920 

0.5 

9 

9 

dp/dT (bar K _1 ) 

-137 

-135 

-134 

5 

10 

10 

Orthobaric densities and temperature of maximun c 

iensity TMD 

TMD/K 

277 

278 

277 

2.5 

10 

10 

P298K/g Cm -3 

0.997 

0.993 

0.99628 

0.5 

9 

10 

P400k/ g cm -3 

0.9375 

0.93 

0.9368 

0.5 

8 

10 

p450K/g c.m “ 3 

0.8903 

0.879 

0.8842 

0.5 

7 

9 

Isothermal compressibility / 10 6 bar 

b 





Kjr [1 bar; 298 K] 

45.3 

46 

45.8 

5 

10 

10 

k t [1 bar;360 K] 

47 

50.9 

49.1 

5 

8 

9 

Gas properties I 

p„[350 K] (bar) 

0.417 

0.13 

0.026 

5 

0 

0 

p„[450 K] (bar) 

9.32 

4.46 

2.64 

5 

0 

0 

B2[450 K] (cm 3 mol - 1 ) 

-238 

-476 

-438 

5 

0 

0 

Heat capacity at constant pressure/cal mol *K 1 1 

C p [liq 298 K; 1 bar] 

18 

21.1 

19.1 

5 

7 

9 

Cp[ice 250 K; 1 bar] 

8.3 

14 

11.9 

5 

0 

1 

Static dielectric constant I 

e [liq; 298 K] 

78.5 

58 

78.3 

5 

5 

10 

e [ 1; 240 K] 

107 

53 

63 

5 

0 

2 

Ratio 

1.36 

0.91 

0.80 

5 

3 

2 

T m -TMD-T c . ratios 1 

T m [I*]/T c 

0.422 

0.394 

0.355 

5 

9 

7 

TMD/T C 

0.428 

0.434 

0.410 

5 

10 

9 

TMD-T m (K) 

3.85 

26 

37 

5 

6 

5 

Densities of ice polymorphs/g cm 3 

p[I/, 250 K; 1 bar] 

0.92 

0.921 

0.919 

0.5 

10 

10 

p[II 123 K; 1 bar] 

1.19 

1.199 

1.196 

0.5 

8 

9 

p[V 223 K; 5.3 kbar] 

1.283 

1.272 

1.275 

0.5 

8 

9 

p[VI 225 K; 11 kbar] 

1.373 

1.38 

1.377 

0.5 

9 

9 

EOS high pressure | 

p [373 K; 10 kbar] 

1.201 

1.204 

1.202 

0.5 

10 

10 

p [373 K; 20 kbar] 

1.322 

1.321 

1.318 

0.5 

10 

9 

Self-diffusion coefhcient/cm 2 s 1 

111 D 278 K 

-11.24 

-11.27 

-11.3657 

0.5 

9 

8 

In D 2 9sa: 

- 10.68 

-10.79 

-10.77 

0.5 

8 

8 

111 D 318 A: 

-10.24 

-10.39 

-10.3 

0.5 

7 

9 

E a k j mol -1 

18.4 

16.2 

16.3 

5 

8 

8 

Shear viscosity / mPa s 1 

77 [1 bar; 298 K] 

0.896 

0.855 

0.863 

5 

9 

9 

rj [1 bar; 373 K] 

0.284 

0.289 

0.307 

5 

10 

8 

Orientational relaxation time / ps I 

Tf" [1 bar; 298 K] 

2.36 

2.3 

2.31 

5 

9 

10 j 

Structure j 

Z 2 (F(Q)) 

0 

8.5 

8.5 

5 

8 

8 

X 2 (overall) 

0 

14.8 

14.8 

5 

7 

7 

Phase diagram 

8 

8 

Overall score (out of 10) 

7.13 

7.31 
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Figure 6: Dielectric constant versus pressure for T — 240,260,280,300,320 K for the TIP4P/£ 
model (circles) and for T = 300 K for experiments 31 (black squares). 


Conclusions 


In this paper we have computed explicitly the behavior of the density and diffusion anoma¬ 
lies for the TIP4P/fi model showing that it gives similar results when compared with the 
TIP4P/2005 model. 

Then, the behavior of the dielectric constant with the temperature and pressure was 
analyzed and compared with experimental results showing better agreement than the values 
obtained by other non-polarizable models. 

Finally, the set of properties proposed by Vega el al. was computed. Our results that 
the performance of the TIP4P/£ model is similar to the performance of the TIP4P/2005 
model 14 with the exception of the dielectric constant for which the TIP4P/fi shows a better 
agreement with the experimental results. We hope that this model might be suitable for 
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studying mixtures and confined systems where the dielectric constant play an important 
role. 
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